Exploring the energy landscape of XY models 



Rachele Nerattini, 1 ' 2 - * Michael Kastner, 3,4, * Dhagash Mehta, 5 '* and Lapo Casetti 1 ' 2,5 

1 Dipartvmento di Fisica e Astronomia and Centro per lo Studio delle Dinamiche Complesse (CSDC), 
Universitd di Firenze, via G. Sansone 1, 1-50019 Sesto Fiorentino (FI), Italy 
2 Istituto Nazionale di Fisica Nucleate (INFN), Sezione di Firenze, 
via G. Sansone 1, 1-50019 Sesto Fiorentino (FI), Italy 
3 National Institute for Theoretical Physics (NITheP), Stellenbosch 7600, South Africa 
4 Institute of Theoretical Physics, University of Stellenbosch, Stellenbosch 7600, South Africa 
5 Department of Physics, Syracuse University, Syracuse, New York 13244, USA 
(Dated: November 21, 2012) 

We investigate the energy landscape of two- and three-dimensional XY models with nearest- 
neighbor interactions by analytically constructing several classes of stationary points of the Hamil- 
tonian. These classes are analyzed, in particular with respect to possible signatures of the thermo- 
dynamic phase transitions of the models. We find that, even after explicitly breaking the global 
0(2) symmetry of the XY spins, an exponentially large class of stationary points are singular and 
occur in continuous one-parameter families. This property may complicate the use of theoretical 
tools developed for the investigation of phase transitions based on stationary points of the energy 
landscape, and we discuss strategies to avoid these difficulties. 
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I. INTRODUCTION 

The stationary points of a function of many variables 
/ : M i— > M. are the points p s G M on a manifold M such 
that V/(p 8 ) = 0. Stationary points play an important 
role for quite a few methods in theoretical physics, as 
knowledge about these points can be used to infer phys- 
ical properties of the system under investigation. When 
the function / is the energy of a many-body system and 
the manifold M is the phase space or the configuration 
space of the system, these methods are referred to as 
energy landscape methods [1]. Examples of applications 
include clusters [1], disordered systems and glasses [2, 3], 
biomolcculcs, and protein folding [4]. Based on knowl- 
edge about the stationary points of the energy function, 
landscape methods can be applied to estimate dynamic 
as well as static properties. In most applications, like 
for example Stillinger and Weber's thermodynamic for- 
malism [5, 6] and other 'superposition approaches' [1, 7] 
for the study of equilibrium properties, only minima of 
the energy landscape are taken into account. In some 
later work, also first-order saddles (sec e.g. Ref. [8]) and 
stationary points of arbitrary index [9] have been consid- 
ered, for instance to characterize glassy behavior [10, 11]. 

The potential relevance of energy landscape properties 
for equilibrium phase transitions was suggested after it 
was realized that stationary points of the Hamiltonian 
are related to topology changes of the phase space ac- 
cessible to the system. It was conjectured that some of 
these topology changes, and therefore some of the sta- 
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tionary points, are at the origin of thermodynamic phase 
transitions [12-16]; quite some research activity followed 
[17-38], some focused on specific models, others trying 
to shed light on the general mechanisms (see [39, 40] for 
reviews) . 

Although equilibrium phase transitions in systems 
with non-fluctuating particle numbers have been mainly 
studied within the canonical ensemble, the connection 
between stationary points of the Hamiltonian and equi- 
librium statistical properties is more transparent in a mi- 
crocanonical setting [41] . This can be understood by ob- 
serving that, for a system with ./V degrees of freedom, the 
entropy density is defined as [42] 



(1) 



where e — E/N is the energy density and uj is the density 
of states. For a system described by continuous variables, 
uj can be written as 



8{U - Ne) dT 



(2) 



where Y denotes the phase space and dY its volume mea- 
sure, S e is the hypersurface of constant energy E = Ne, 
and d£ stands for the N — 1-dimensional Hausdorff mea- 
sure. The rightmost integral stems from a co-area for- 
mula [43]. At a stationary point p s , the gradient VH(p s ) 
vanishes by definition and the integrand diverges, while 
at the same time the measure dE shrinks such that ui 
in general remains finite for finite systems. Indeed, a 
more refined analysis [44] shows that, although the inte- 
gral on the right-hand side of (2) remains finite in the 
vicinity of a stationary point, the density of states will 
be nonanalytic at stationary values e s := H(p s )/N of the 
energy density for any finite N [45]. The microcanonical 
nonanalyticities appearing at finite N are found to be 
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in correspondence with stationary configurations; how- 
ever, the 'strength' of such nonanalyticities generically 
decreases linearly with N, i.e., the first k derivatives of 
the entropy are continuous, where k is O(N) [44, 46]. The 
usual thermodynamic quantities, like equations of state, 
are given by low-order derivatives of the entropy, and the 
observation of nonanalyticities of order O(N) from noisy 
data is therefore restricted to very small system sizes N. 
Taking this to its logical conclusion, we expect the order 
of the nonanalyticities to diverge in the thermodynamic 
limit, leading to a vanishing effect of stationary points 
and smooth results for the thermodynamic functions. 

Superficially, this seems to challenge the conjecture 
that phase transitions stem from stationary points of 
the energy, as it suggests that finite- N nonanalyticities 
due to stationary points arc unrelated to thermodynamic 
phase transitions [47] . And indeed, for several model sys- 
tems, phase transitions have been found to occur at ener- 
gies at which no stationary points of the Hamiltonian are 
present [18, 21, 22, 24, 26, 29, 48, 49]. On the other hand, 
a substantial amount of evidence (in the form of model 
calculations) has accumulated in favor of the conjecture 
that stationary points often do play a relevant role for 
the emergence of phase transitions, and the presence of 
a transition reflects prominently in properties of the sta- 
tionary points [50]. This evidence comes mostly from 
exactly solvable systems (often with mean-field interac- 
tions) where the connection between stationary points 
and thermodynamic phase transitions has been shown 
explicitly [13, 15-17, 23, 51]. Subsequently, and of di- 
rect relevance to the present work, a possible scenario 
of how certain finite- N singularities may survive in the 
thermodynamic limit has been proposed in [52, 53]: Only 
those singularities related to asymptotically flat station- 
ary points may survive in the thermodynamic limit and 
induce a thermodynamic phase transition. Asymptoti- 
cally flat here refers to stationary points whose deter- 
minant of the Hessian matrix of the potential energy V 
vanishes in the thermodynamic limit. More precisely, for 
a sequence of potential energy functions Vn : Mjv h> R 
on configuration space Mat, consider a sequence of sta- 
tionary points {q s N G Mn} N£K such that VViv(^) = 0. 
If the limit 

v c = lim V N {q s N )/N (3) 

N— >oo 

exists and 

lim |detHess Vjv (q s N )\ 1/N - 0, (4) 

N— >oc 

then the stationary points may induce a singularity in 
the configurational entropy density s(v) at v — v c in the 
thermodynamic limit. This indeed was verified to happen 
in the exactly solvable mean-field models where the con- 
nection between stationary points of the Hamiltonian and 
phase transitions had been previously demonstrated [53] , 
and also in a non-solvable toy model of a self-gravitating 
particles with a phase transition between a homogeneous 
and a collapsed phase [31]. Note that several of the above 



results, including the determinant criterion (4), are de- 
rived under certain gcnericity assumptions on the Hamil- 
tonian or potential energy function. In particular, all 
stationary points are required to be isolated, and we will 
come back to this condition in Sec. VI. 

The main purpose of the present paper is to go beyond 
one-dimensional or mean-field systems by performing an 
analysis of stationary points and their Hessian determi- 
nants for classical XY spin models with nearest-neighbor 
ferromagnetic interactions. In particular, we consider the 
XY model on a two-dimensional square lattice and on a 
three-dimensional cubic lattice. These models arc among 
the simplest lattice spin models with short-range inter- 
actions amenable of an energy landscape approach based 
on stationary points of the Hamiltonian (in the even sim- 
pler Ising model such an analysis is impossible due to 
the discrete character of the spin variables) . For the one- 
dimensional XY model [34] as well as for the mean-field 
(fully connected) XY model [16, 53] the stationary points 
and their relation to phase transitions have already been 
studied in earlier work. In both cases a fully analytical 
treatment of all the stationary points was feasible, as well 
as an analysis of flat stationary points in the sense of Eq. 
(4). The results showed a clear signature in stationary- 
point properties of the presence of a finite-temperature 
phase transition in the mean-field XY model; and no sig- 
nature of a transition in the one-dimensional XY model, 
in agreement with the known thermodynamic behavior 
of these models. 

In the more interesting case of short-range interacting 
XY models in dimensions d 2, it seems unlikely that 
finding all stationary points of the Hamiltonian is fea- 
sible, at least for lattices of reasonable sizes (for small 
lattices, the stationary point analysis has been done in 
Refs. [54-56] using a homotopy continuation method). 
It is, however, possible to analytically construct certain 
classes of stationary points of the Hamiltonian, and this 
will be the main topic of the present article. Two of 
these classes of stationary points are large in the sense 
that they contain an exponentially (in the lattice size 
N) growing number of points. Exponential growth with 
N is also expected for the overall number of stationary 
points [5, 57, 58] and in this sense we may at least hope 
that the classes of stationary points we have constructed 
are in some way relevant for the thermodynamics of the 
model. We evaluated the Hessian determinant at some 
of these stationary points and tried to identify candidate 
sequences of stationary points satisfying the criterion (4). 

After introducing the nearest-neighbor XY model in 
Sec. II and deriving expressions for stationary points and 
Hessian matrices in Sec. Ill, a first class of what we call 
Ising stationary configurations is constructed in Sec. IV: 
these are stationary points of the XY model where each 
pair of neighboring spins is either aligned or anti-aligned. 
Although these configurations are easy to construct, we 
have to resort to sampling strategies to select members of 
this exponentially large class, and evaluate the Hessian 
determinant numerically Earlier work [59, 60] suggests 
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that these Ising stationary configurations might be par- 
ticularly relevant for the phase transition of the model. 
However, unlike in the mean-field case, no signature of 
the phase transitions can be discerned from the data we 
obtained for this class of stationary points for the two- 
dimensional and three-dimensional models. In Sec. V wc 
construct a class of stationary points which have the char- 
acter of spin waves. Since their number is small (subcx- 
ponential), they are not expected to significantly influ- 
ence the thermodynamic behavior of the model. In Sec. 
VI, we construct an exponentially large class of particu- 
larly interesting stationary points whose Hessian deter- 
minant is zero. Moreover, we are able to prove that, even 
after explicitly breaking the global 0(2) symmetry of the 
XY model, these stationary points are not isolated, but 
occur in continuous families. This finding has interest- 
ing consequences which will be discussed at the end of 
Sec. VI. In Sec. VII we investigate how the presence of 
inhomogeneous external magnetic fields may destroy the 
continuous families and lead to isolated, nonsingular sta- 
tionary points. A summary of the results and concluding 
remarks are presented in Sec. VIII. 

II. NEAREST-NEIGHBOR XY MODELS 

A paradigmatic class of models for the study of mag- 
netic phase transitions (the prototype of all continuous 
phase transitions) are classical 0(n) spin lattice mod- 
els. We consider d-dimensional hypercubic lattices with 
periodic boundary conditions. To each lattice site i an 
n-component classical spin vector Si = (Sj,...,S^) of 
unit length is assigned. The energy of the model is char- 
acterized by the Hamiltonian 

n 

ft(«) = - J J2 Si ■ Sj = - J e E s t s h ( 5 ) 

(ij) a=1 

where the angular brackets denote a sum over all pairs of 
nearest-neighboring lattice sites. The exchange coupling 
J will be assumed to be positive, resulting in ferromag- 
netic interactions, and without loss of generality we set 
J = 1 in the following. The Hamiltonian (5) is globally 
invariant under the 0(n) group. For the case n = 1 the 
symmetry group 0(1) = Z 2 is a discrete group and the 
Hamiltonian (5) becomes the Ising Hamiltonian 

ft«=-5>^ (6) 

(i,3) 

where cr, G { — 1, +1} Vi. In all other cases n ^ 2 the 0(n) 
group is continuous. As special cases, the 0(n) models 
include the XY model for n — 2 and the Heiscnbcrg 
model for n = 3. We will consider in the following the 
XY model on a two-dimensional square lattice and on 
a three-dimensional cubic lattice, hoping that this choice 
yields the simplest possible 0(n) models amenable to the 
energy landscape analysis we have in mind. The XY 



spins live on the unit circle S 1 , so that the components 
of the zth spin can be parametrized by a single angular 
variable $j G [0, 2tt), 

is} =cos^, 

We can thus conveniently write the Hamiltonian of the 
XY model as 

1 N 

^^E E ( 8 ) 

»=i 

where Af(i) denotes the set of nearest neighbors of lattice 
site i. The energy density e = VS 2 ^ /N takes values in the 
interval [— d, d] where d is the lattice dimension. 

In two dimensions, the XY model exhibits a Bere- 
zinskij-Kosterlitz-Thoulcss (BKT) phase transition [61, 
62] with no long-range order, while in three dimensions it 
undergoes a ferromagnetic transition in the same univer- 
sality class as that of the superfluid transition [63] . With 
our choices for the lattices and the interaction constants, 
the BKT transition in the two-dimensional case occurs 
at a critical temperature T c 2d « 0.8929(1) correspond- 
ing to a critical energy density e 2d — (H^) (T c ) /N w 
— 1.4457(4) [64, 65], while the critical temperature of the 
ferromagnetic transition in three dimensions is T c 3d « 
2.20167(9), corresponding to a critical energy density 
£ 3d K _ .99184(6) [66]. 

III. STATIONARY POINTS AND HESSIAN OF 
THE ENERGY 

To conduct an energy landscape analysis as outlined in 
the Introduction, we have to find the stationary points of 
the energy, i.e., solutions $ s satisfying the vector equation 
VH (2) (i? s ) = 0. Using (8), the fcth component of this 
equation can be written as 

E sin(0 fc -0 3 -)=O. (9) 

The 0(2) invariance of the Hamiltonian (8) implies that 
the solutions of (9) are not isolated points in configu- 
ration space, but occur in continuous curves: Given a 
stationary point $ s = , . . . , $ S N ), the points (#f + 
a, . . . , i!) s N + a) are also stationary points for arbitrary 
a G R. Several of the theoretical tools and results men- 
tioned in the Introduction, and in particular the Hessian 
determinant criterion (4), require energy functions with 
only isolated stationary points. It is therefore necessary 
to explicitely break the global 0(2) symmetry of the XY 
model. Here we choose to fix one spin, e.g. An = 0, but 
there are other possibilities how the symmetry breaking 
can be achieved, and we will come back to this point in 
Sec. VII. In the large- N limit, the only effect of this global 
phase fixing is to dictate the direction of the spontaneous 
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symmetry breaking, but otherwise thermodynamic quan- 
tities remain unaffected. 

In order to apply the criterion (4), we have to evaluate, 
at the stationary points, the determinant of the Hessian 
matrix of the Hamiltonian [which for the XY model (8) 
coincides with the potential energy] . The elements of the 
Hessian matrix are defined as 



d 2 n 



(10) 



The constraint = makes the Hessian an (N — 1) X 
(N — 1) matrix and, for the XY Hamiltonian (8), its 
diagonal elements are given by 



H 



kk 



y] cos 

jeAf(fc) 



(11) 



while the off-diagonal elements are 



Hui = 



- cos (i? fc - forZeTV(fc), 
else, 



(12) 



for k,l = 1,...,N- 1. 



IV. ISING STATIONARY CONFIGURATIONS 

Finding all stationary points of the Hamiltonian (8) is 
unlikely to be feasible for large lattices. This section and 
the subsequent Sees. V and VI will be concerned with 
the construction of special classes of stationary points 
and with the evaluation of their Hessian determinant. 

Inspection of the stationary point conditions (9) re- 
veals that any configuration where = {0, it} Mi is a 
stationary point, as in this case each term of the sum 
on the left-hand side of (9) vanishes separately. In the 
notation of (7) such stationary points can be written as 



S} — 



0. 



(13) 



where Ui E {— 1,+1}. Therefore, as already discussed 
in [59], each such stationary point -d s of the XY Hamil- 
tonian (5) corresponds to a configuration of the Ising 
model (6) defined on the same lattice. Moreover, the 
corresponding stationary values TL^(i3 s ) of these 'Ising 
stationary configurations' are just the energy levels of the 
Ising Hamiltonian (6). 

As we will see in later sections, the Ising stationary 
configurations are not the only stationary points of the 
XY model, but they constitute a large class in the sense 
that their number grows like 2 N with the system size 
N. But, as already mentioned, also the overall num- 
ber of stationary points of a generic function of N vari- 
ables is expected to scale exponentially with N, a fact 
that may be seen as a hint towards the relevance of the 
Ising stationary configurations. Moreover, as discussed 



in [59, 60], the critical energy densities of O(n) models 
are remarkably close to those of the corresponding Ising 
model, another indication pointing towards the relevance 
of the Ising configurations for the phase transition of the 
XY model. 

Evaluated at the Ising stationary configurations, the 
Hessian matrix elements (11) and (12) can be written as 



Hkk = Cfe 



jeAf(k) 



and 



Hm = 



— o-fccrj for Z eJV(fc), 
else, 



(14) 



(15) 



for k, I — 1, . . . ,N — 1 and k ^ I. In the mean-field 
XY model [16], and also in other mean-field or one- 
dimensional models [31, 34, 53], the energy and the Hes- 
sian determinant of Ising stationary configurations de- 
pend on only a single collective variable, thus allowing an 
analytical search of stationary points satisfying (4). Un- 
fortunately, for the two- and three-dimensional nearest- 
neighbor models this is not the case and we have to resort 
to numerical methods. We computed the determinant of 
the Hessian of the Hamiltonian on a numerically obtained 
sample of the Ising stationary configurations. The sample 
was obtained by standard Metropolis Monte Carlo sim- 
ulations of the two-dimensional and three-dimensional 
Ising models, exploiting the above mentioned one-to-one 
relation between configurations of the Ising model and 
Ising stationary configurations of the XY models. 



A. Two-dimensional XY model 

We considered L x L square lattices of side lengths 
L = 16, 24, 32 and 64, so that the number of degrees of 
freedom ranges from N = L 2 = 256 to 4096. Compared 
to those typically considered in simulations nowadays, 
these are not very big lattices, and indeed obtaining the 
sample was easy and fast. The practical limit on the num- 
ber of degrees of freedom was set by the time-consuming 
calculation of the Hessian determinant for each configu- 
ration of the sample. Although in principle Ising configu- 
rations occur over the entire range [—2, 2] of accessible en- 
ergy densities, only configurations with negative energy 
were sampled in the Monte Carlo runs. This is a con- 
sequence of using canonical simulations at positive sim- 
ulation temperature so that, for sufficiently large lattice 
sizes, the Boltzmann weight narrowly focuses the sam- 
pled distribution on a range of negative energies. How- 
ever, by using also negative temperatures we would have 
obtained symmetric results with respect to zero energy, 
without adding any relevant information. For each lat- 
tice of side lengths L = 16, 24 and 32, we considered a 
total sample of 250000 configurations. For L = 64 we 
considered only 48000 configurations, the Hessian deter- 
minant being quite heavy to compute. Results for the 
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FIG. 1. (Color online) Rescaled Hessian determinant D of 
Ising stationary configurations for the two-dimensional XY 
model, plotted as a function of the energy density e. Data 
symbols correspond to lattices of side lengths L ranging from 
16 to 64 (see the legend in the plot for the color/grayscale 
code). The critical energy density e\ d « -1.446 of the BKT 
transition is marked by a vertical dashed line. The solid lines 
are the values calculated for the polygonal configurations in 
the large- N limit according to (24). 
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FIG. 2. (Color online) As in Fig. 1 for the index density i. 
Data for the polygonal configurations are not shown. 



rescaled Hessian determinant 

D = |detHess w( 2) (tf s )\ 1/N 



(16) 



as a function of the energy density are shown in Fig. 1. 

In order to further characterize the sampled station- 
ary points we computed, for the same lattices, the index 
density 



index(?9 s 
N - 1 



(17) 



where the index of a stationary point f? s is the number 
of negative eigenvalues of the Hesse matrix at $ s . The 
results for the index density l as a function of the energy 
density are shown in Fig. 2. 

Two features of the results shown in Figs. 1 and 2 are 
of particular interest. 



(i) As N grows, the rescaled determinant D as well 
as the index density i show a tendency to concen- 
trate onto a single curve, so that, at least for Ising 
stationary configurations, these quantities appear 
to be good thermodynamic observables. Moreover, 
both quantities appear to be monotonic functions 
of the energy density. 

(ii) The Hessian determinant shows no tendency to 
vanish for any value of the energy density. Hence 
there are no indications of the presence of asymp- 
totically flat stationary points, i.e., of the valid- 
ity of Eq. (4) around the transition energy density 
e\ d « -1.446 of the BKT transition. Also the index 
density i(e) does not show any remarkable feature 
close to e1 d . 

Our sample has variable magnetization and, in particu- 
lar for low energies, configurations typically have nonzero 
magnetizations while in the two-dimensional XY model 
the typical magnetization is zero at any energy. In or- 
der to rule out the possibility that this may affect our 
results, we repeated the calculation of the Hessian de- 
terminant on a sample of configurations with vanishing 
magnetization, obtained by Monte Carlo with Kawasaki 
dynamics. The results (not shown) display no apprecia- 
ble differences with respect to Fig. 1. 



B. Three-dimensional XY model 

In the three-dimensional case we proceeded analo- 
gously to the two-dimensional case, considering LxLxL 
lattices of side lengths L = 8, 10 and 12, so that the num- 
ber of degrees of freedom ranged from N = L 3 = 512 to 
1728. For each lattice we considered a total sample of 
57000 configurations. Results for the rescaled Hessian 
determinant D and for the index density t as a function 
of the energy density are shown in Figs. 3 and 4. The 
similarities to Figs. 1 and 2 are striking, so that the con- 
siderations made for the two-dimensional case carry over 
to three dimensions. 



V. POLYGONAL STATIONARY POINTS 

Another class of stationary configurations of the XY 
model that can be easily identified are those for which 
neighboring spins differ by the same angle tp [67], 



0j = ft ± ip Vi = 1, . . . , N, j e Af(i). 



(18) 



Periodic boundary conditions restrict these angles to val- 
ues if = 27rm/L with m € {0, . . . , L — 1}. The stationary 
configurations (18) are analogous to those called 'spin- 
wave' stationary points for the 1-d XY model in [34] and 
to the 'polygonal' stationary points of the self-gravitating 
ring model introduced in [31]. We will adopt the latter 
naming convention. 
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FIG. 3. (Color online) Rescaled Hessian determinant D as 
a function of the energy density e for the three-dimensional 
XY model. Data symbols correspond to lattices of side 
lengths L ranging from 8 to 12 (see the legend in the plot 
for the color/grayscale code). The critical energy density 
ej? d ?a —0.99 of the ferromagnetic transition is marked by a 
vertical dashed line. The solid lines are the values calculated 
for the polygonal configurations in the large- TV limit according 
to (25). 




FIG. 4. (Color online) As in Fig. 3 for the density of index . 
Data for the polygonal configurations are not shown. 



For the two-dimensional XY model, the energy density 
of a polygonal stationary configuration is 



e{(p) = —2 cosy, 



(19) 



and the Hessian determinant of the Hamiltonian has the 
simple form 

Hij (<p) = Ai,j cos (p, (20) 
where A is an (N — 1) X (N — 1) matrix with elements 

4 i£i=j, 



A i:j = ^ -1 if j G Af(i), 



(21) 



else. 

We have analyzed the determinant of A numerically, and 
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FIG. 5. (Color online) Rescaled determinant D = 
| det .Al 1 /^ -1 ' of the matrix A as defined in (21), plotted as a 
function of the inverse system size. The line is obtained from 
a linear least-square fit, and an extrapolation to 1/L 2 = 
yields a = limjv-^oo D ~ 3.21. 



the results shown in Fig. 5 provide strong evidence that, 
asymptotically for large N, the determinant behaves as 



dot A 



„JV-1 



(22) 



with a « 3.21. The rescaled Hessian determinant com- 
puted on these configurations in the thermodynamic limit 
is then given by 



lim \det H((f)\ 

N— >oo 



1/(JV-1) 



and, using (19), we can write 



N 



lim \dctH((p) 



i/(N-l) 



a\ cos ip\ 



a \ i 



(23) 



(24) 



This result is plotted in Fig. 1 along with the data for 
the Ising stationary points. 

For the polygonal stationary points of the three- 
dimensional XY model, the calculation proceeds along 
very similar lines, yielding as a final result 



lim Idet I 

N— >oo 



1/(JV-1) 



(25) 



with b w 5.3. This result is plotted in Fig. 3, along with 
data for the Ising stationary configurations. 



VI. SINGULAR STATIONARY POINTS 

In the Introduction we elaborated, amongst other 
things, on the criterion (4), which relates the occurrence 
of a phase transition in the thermodynamic limit to the 
existence of stationary points whose Hessian determi- 
nant, in this limit, approaches zero in a suitable way. 
This, and other results mentioned, were derived under 
the assumption that the Hamiltonian H of the system 
under consideration is a Morse function, meaning that 
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at any stationary point of H the Hessian determinant is 
nonzero. In principle, this requirement is not particularly 
restrictive: Morse functions are known to be generic, in 
the sense that, even if a given function is not a Morse 
function, it can be made into one by adding a generic 
perturbation (see Chapter 4.4 and also the introductory 
comments in Chapter 3.1 of [68]). Although this prop- 
erty underlines the relevance of Morse function for realis- 
tic (and therefore imperfect) physical systems, the Morse 
property may well be violated in clean, idealized physical 
models. To avoid this problem, one could of course add 
a small generic perturbation to the model, but this will 
typically make analytic calculations much harder, if not 
impossible. 

In the following we prove that, in lattice dimensions 
d = 2 and d = 3, the XY Hamiltonian (8) is not a Morse 
function, but instead has an exponentially (in N) large 
number of singular stationary points. Moreover, the sta- 
tionary energy densities 'H('d s )/N of all these singular 
stationary points become dense on the interval [—d, d] of 
accessible energy densities in the thermodynamic limit. 
The proof is constructive, and for simplicity we restrict 
the presentation to two-dimensional square lattices of size 
L x L with periodic boundary conditions. The three- 
dimensional case is treated in Appendix A. Generaliza- 
tions to higher-dimensional lattices should be possible 
along similar lines, but we did not work this out in de- 
tail. 

For a configuration $ = . . . , to have a vanish- 
ing Hessian determinant, it is sufficient that one row of 
the Hessian matrix given in (11) and (12) has only zero 
entries. This is a local property, as all the nonzero en- 
tries in the fcth row are fully determined by the fcth spin 
and its nearest neighbors. Consider for example a con- 
figuration which, somewhere on the lattice, locally looks 
like 

• I ' 

t <- f (26) 

• I ' 

where arrows f, — J., <— correspond to angle variables 
$i = 0, 7r/2, 7T, 3ir/2. The dots in (26) are place hold- 
ers for arbitrary spin orientations, as their values do not 
matter for the moment. Assigning to the center (left- 
pointing) spin of this configuration the label k, we find 

H kk = cos(±7r/2) - 0, H kl =0Vl (27) 

for the elements of the Hessian matrix. The matrix there- 
fore does not have full rank and its determinant is zero. 

Next, in addition to this local condition guaranteeing 
that the Hessian determinant vanishes, we also have to 
ensure that the overall configuration is a stationary point 
of T-i^ 2 \ This is a global property as, in order for a con- 
figuration to be stationary, the constraint 

J2 frin(<?fc-i?i) = (28) 



has to be satisfied for all lattice sites k. Starting from 
the pattern in (26), it is not too difficult to construct an 
embedding of such patterns into larger lattices while at 
the same time satisfying the stationarity constraints (28). 
For the example of an 8 x 8 square lattice, the following 
class of configurations does the job, 
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X 
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t 
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X 
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-> 
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X 


X 


X 


X 


I 



(29) 



The lattice sites marked with gray ^-arrows can be filled 
with an arbitrary 'Ising-typc'-pattcrn of t and 4- arrows. 
Independently of the precise pattern of these up- and 
down-pointing arrows, the resulting configuration will al- 
ways be stationary. In this way, we have obtained a class 
of stationary points of the Hamiltonian Ti^ with vanish- 
ing Hessian determinant, and the scheme works in just 
the same way for larger lattice sizes. 

This class of singular stationary points i? s is ample 
enough to allow us to adjust the energy density H(d s )/N 
almost freely: By choosing an appropriate Ising-type pat- 
tern of | and I for the gray ^-arrows in (29), the en- 
ergy of the configuration is varied. Since the number 
of gray ^-arrows in such a configuration scales as L 2 , 
their contribution to the overall energy will, in the large- 
ly limit, dominate over the fixed (black) arrows in (29) 
whose number increases only linearly in L. As a result, 
the corresponding stationary energy densities H(i9 s )/iV 
are dominated by the Ising-type pattern chosen for the 
gray ^-arrows and, like the Ising energy densities, become 
dense on the interval [—2, 2] of accessible energy densities 
in the thermodynamic limit. 

Singular stationary points come in two flavors: They 
can either be isolated stationary points, like at the min- 
imum x s = of the quartic f\(x) = x 4 . Or they 
can form continuous families of non-isolated stationary 
points, like for the Mexican hat potential f2{x,y) = 
(x 2 + y 2 ) 2 — (x 2 + y 2 ) where the points on the circle 
x 2 + y 2 = 1/2 form a continuous curve of minima of 
f-2 ■ Our singular stationary points of the two-dimensional 
XY Hamiltonian fall into the latter category. This can 
be seen by starting from a configuration like the one de- 
picted in (29) and then simultaneously rotating by some 
arbitrary angle a all the — > and <— spins situated on the 
diagonal. It is easily checked that the resulting configu- 
ration still satisfies the stationarity condition (28). This 
proves that the singular stationary points we have con- 
structed are not isolated, but occur in continuous one- 
parameter families, parametrized by the angle a. Sim- 
ilarly, one can create two- and more-parameter families 
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by generalizing (29) to contain more than one diagonal 
pattern, 



t 


t 


t 




t 


t 


t 




X 


X 


-> 




X 


; 
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t 




t 


$ 


-> 
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1 



(30) 



In this configuration the — > and ^— spins situated on the 
main diagonal can be simultaneously rotated by some 
angle a, and those on the other diagonal (modulo peri- 
odic boundary conditions) by an independent angle /3, 
resulting in a continuous two-parameter family of sta- 
tionary points. The generalization to more parameters 
is straightforward, provided the lattice sizes are chosen 
large enough. 

Note that this occurrence of continuous families of non- 
isolated stationary points is not due to the global 0(2) 
invariance of the XY Hamiltonian: This global symme- 
try is a trivial effect that we have taken care of by fixing 
one angle variable, = 0, as discussed in Sec. III. Prom 
the examples (29) and (30), however, we have learned 
that this global phase fixing is not sufficient to ensure 
that the XY Hamiltonian is a Morse function with only 
isolated stationary points. The problem seems to be that 
certain spin environments, like the pattern 



t 



■ • • • t ■ 

• • • X ■ X 

■ • t • t 

■ X ■ X 

t ■ t • 

X ■ X • ■ 

t • t ■ • • 



(31) 



X 



in (29), can build a 'cage' around a lattice region such 
that the overall phase of the enclosed region [the diago- 
nal in the case of (29)] is shielded from the rest of the 
configuration. As a consequence, breaking of the global 
0(2) invariance of the XY Hamiltonian by locally fix- 
ing $jv = is not sufficient. Another way to eliminate 
the global 0(2) invariance is to use antiperiodic bound- 
ary conditions in all the d-directions, as proposed in Ref. 
[69]. However, we have verified numerically that even 
using antiperiodic boundary conditions, isolated singular 
solutions as well as continuous one- and more-parameter 
families of solutions exist. 



VII. SYMMETRY BREAKING FIELDS 

The observation of local versus global properties also 
suggests how the problem of non-isolated, singular sta- 
tionary points might be solved: As mentioned at the be- 
ginning of Sec. VI, perturbations like 



N 



N 



U(2) = -\J2 12 cos (0< -<?,•)- ^ M< (32) 



i=l 



and maybe also 



JV 



n 



(2) 



\ll 12 cos ^ 



JV 

12 

i=l 
>N 



hi cos di, (33) 



for generic values of (hi , ■ ■ ■ , hjv) & Mr, should ensure 
that the Hamiltonian has only isolated and nondegener- 
ate stationary points, but other forms of perturbations 
might do the job as well. For 3x3 lattices we have 
checked numerically that, up to numerical accuracy, the 
perturbations in (32) and (33) indeed destroy all singular 
stationary points of Ti^: Firstly, we used the numeri- 
cal polynomial homotopy continuation method (Bertini 
software package [70]) which finds all the solutions of a 
system of multivariate polynomial equations, including 
isolated singular solutions [71]. This method has been 
recently used to study the potential energy landscape in 
various areas of physics [35, 54-56, 72, 73]. We studied 
at least 10 generic sets of (hi,..., ft. at) £ ~M N for both 
types of perturbations and verified that no isolated sin- 
gular solution occur for these perturbed systems. We 
then used an extension of the numerical polynomial ho- 
motopy continuation method, called numerical algebraic 
geometry [74, 75] , which can find solution curves of a sys- 
tem of polynomial equations, combined with the method 
described in Ref. [76], and concluded that there is no 
continuous solution curve for any of the systems in the 
presence of a generic perturbation. 

From a physical point of view, the cosine-perturbed 
Hamiltonian (33) appears particularly appealing as it has 
the form of a spatially inhomogeneous magnetic field in 
x-direction acting on the spins. Interestingly, this spe- 
cific choice of the perturbation leaves the Ising stationary 
configurations of Sec. IV unaffected: Every Ising config- 
uration ($f, . . . , $^f) with € {0, 7r} is also a stationary 
point of the perturbed Hamiltonian (33) for arbitrary 
perturbation fields hi . Mathematically, this is due to the 
fact that the Taylor expansion of the perturbation around 
an Ising stationary configuration 



JV 

hi cos 

i=i 



JV 



hi cos0J + 0(0i -^) S 



(34) 

has vanishing linear contributions, thus leaving these sta- 
tionary points unaffected. It is unclear to the authors 
whether there is any physical significance to this obser- 
vation. This notwithstanding, this property can be used 
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to check if all the singular solutions of are indeed 
destroyed by the perturbation (33) also in lattices larger 
than 3x3. For simplicity in the following we will restrict 
ourselves to two-dimensional square lattices, but we have 
checked that the conclusions remain valid in three dimen- 
sions. 

In order to study the effect of the perturbation (33) on 
singular configurations, we want to construct a sample 
of such configurations, spread over a range of energies 
similar to the nonsingular ones in Fig. 1. One possible 
strategy to do so is to take the nonsingular sample as a 
starting point and transform each of the configurations 
into a singular one by imprinting the mask 



■■••It 

• ■ -tit 

• ' I t I 
"tit 
I f I ■ 

t | t • ' 
I t I • • • 
I t • ■ • • 



(35) 



t 



(or a similar one for other lattice sizes), i.e., by rotating 
all spins of the configuration into the orientation indi- 
cated in (35), while leaving unchanged all sites indicated 
by dots. The configuration in (35) is similar to the one 
in (29), only that the spins on the diagonal are rotated 
by 7r/2. Such a configuration, as explained in Sec. VI, is 
also singular, and it preserves the Ising-type character of 
the configuration. Imprinting the mask (35) causes only 
a subextensivc change of energy, and the distribution in 
energy density of the stationary points will therefore be 
similar to the distribution of the original (nonsingular) 
sample. Switching on the perturbation fields hi in (33) 
should turn all singular solutions into regular ones, and 
it is this effect we want to study. 

We performed the above analysis on 25000 configura- 
tions for a square lattice of side length L = 24. The fields 
hi were chosen randomly in the ranges ri = [—0.5, 0.5] 
and r 2 = [— I0 -7 , 10~ 7 ], to test the dependency of the 
reduced determinant on the strength of the fields hi . Re- 
sults are shown in Fig. 6. This analysis, like the one 
conducted for the 3x3 lattice by numerical homotopy 
continuation, confirms that generic perturbations as in 
(33) transform singular solutions of into nonsingu- 
lar ones. Remarkably, the effect of the perturbations hi 
on the rescaled determinant D is rather drastic: Already 
for tiny perturbations in the range r 2 = [— 10 -7 , 10 -7 ], 
D is far away from zero and very close to the values of 
the original (nonsingular) Ising stationary configurations. 
This finding can be explained by the fact that, according 
to the scheme in Sec. VI, we constructed one-parameter 
families of singular solutions. Accordingly, the Hessian 
determinant is expected to have a single vanishing eigen- 
value. Switching on a perturbation affects the zero eigen- 




FIG. 6. (Color online) Rescaled Hessian determinant D as 
a function of the energy density for the two-dimensional XY 
model with L = 24 and cosine perturbation terms (33). See 
the legend in the plot for the color /grayscale code: "origi- 
nal" stands for the original Ising configurations as in Fig. 1; 
"modified" stands for the singular Ising configurations (built 
from the original ones as described in the text), either with 
unperturbed Hamiltonian, so that D — 0, or with a perturbed 
Hamiltonian with the fields hi chosen randomly in the range 
n = [-0.5, 0.5] or r-2 = [-10~ 7 , 10" 7 ]. 



value by making it nonzero and of order h, while all other 
eigenvalues remain constant (nonzero) to leading order. 
We therefore have 



D(h) = \ch\ 



l/N 



N 



n ^ 

k=2 



l/N 



(36) 



where the eigenvalues are independent of h to lead- 
ing order for all k 2. No matter how small ch is, 
|c/i| 1/Ar will always be close to 1 for A ^ 1. The remain- 
ing (A — l)-fold product in (36) will generically yield 
a value close to the 'thermodynamic' one observed for 
generic Ising stationary configurations as shown in Fig. 
1, even for tiny perturbations h. Intuitively, we would 
expect a stationary point with an extensive number of 
vanishing eigenvalues to be more relevant for the sys- 
tem's thermodynamic properties, while those with a few 
such eigendirections should not play a major role. But 
this is speculation going beyond what the determinant 
criterion (4) claims and needs further examination. 

In summary, we find that a generic perturbation as in 
(32) or (33) successfully destroys all singular stationary 
points. Moreover, the rescaled Hessian determinant D 
is rather insensitive to the actual strength of the per- 
turbation. Similar behavior is observed for the three- 
dimensional XY model, but the results are not shown 
here. 



VIII. CONCLUSIONS 

We have explored the energy landscape of the XY 
model with nearest-neighbor interactions on the two- 
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dimensional square lattice and the three-dimensional cu- 
bic lattice. In particular, we have constructed certain 
classes of stationary points of the Hamiltonian (8). One 
of these classes consists of Ising stationary configura- 
tions (13), and their number is 2 N for a given lattice 
size N. While analytic expressions for all these expo- 
nentially many stationary points are readily obtained, 
an analysis of their properties is a much harder task. 
We resorted to Monte Carlo techniques for generating 
samples of Ising stationary configurations and then nu- 
merically calculated properties like the index i and the 
rescaled Hessian determinant D of these points. The re- 
sults, summarized in Figs. 1-4, indicate that D and I are 
good thermodynamic observables in the sense that, with 
increasing lattice size N, the data points concentrate on 
a line in these plots and appear to be functions of the 
energy density e alone. 

The original motivation for undertaking this energy 
landscape study was to test whether the criterion (4), 
based on the Hessian determinant at stationary points of 
the Hamiltonian, reveals a signature of the phase tran- 
sition of the XY model in two or three dimensions. In 
this respect, our results are not conclusive. The data 
for the rescaled Hessian determinant D, shown in Figs. 1 
and 3, are clearly bounded away from zero for all values 
of the energy density e, and therefore do not signal the 
presence of a phase transition according to the criterion 
(4). As far as the validity of (4) is concerned, however, 
this finding has little to say. It rather reveals the limita- 
tions of the numerical method we have been using: The 
Monte Carlo technique we have been using to generate 
a sample of Ising stationary configurations uses impor- 
tance sampling with respect to the energy, resulting in 
a reasonably uniform distribution of data points on the 
energy axis in Figs. 1-4. But for a given energy density e, 
stationary points are selected unbiased, resulting in the 
above mentioned behavior as 'good thermodynamic ob- 
servables'. This implies, however, that stationary points 
with vanishing (or at least small) rescaled Hessian de- 
terminant D are found by this sampling technique only 
in case that D = is the most probable value at some 
energy density e (see [77] for a numerical study of the 
nearest-neighbor c/> 4 model on the square lattice reaching 
similar conclusions). According to our data, this is not 
the case. 

Indeed, and rather surprisingly, we were able to show 
that singular stationary points, i.e., stationary points 
with D = 0, do exist and are even in abundance: As 
proved in Sec. VI, even after breaking the global 0(2) in- 
variance of the XY model by fixing one spin, an exponen- 
tially (in N) large number of singular stationary points 
exists, densely covering the accessible range [—d, d] of en- 
ergy densities in the large- TV limit. Moreover, these sin- 
gular stationary points are non-isolated, i.e., they come 
in continuous families parametrized by one or several an- 
gular variables. But despite their ubiquitous presence 
and abundance, our Monte Carlo scheme failed to detect 
these points, as the value D = of their rescaled Hes- 



sian determinant is not the most probable one at any 
given s. It must be noted that this is not a limitation 
of the specific Monte Carlo technique we used here: it is 
expected to be a generic property of unbiased numerical 
sampling schemes. For instance, also a search of sta- 
tionary points by means of a modified Newton-Raphson 
method analogous to that used in Rcf. [77] did not reveal 
any tendency of D to vanish close to the phase transition 
energies but did not find any singular solutions either 
(we have not shown these data in the paper). Hence 
our results suggest that from a practical point of view 
a purely numerical approach to the criterion (4) is not 
very useful unless a numerical sampling scheme able to 
efficiently detect stationary configurations with zero — or 
at least small — determinant is devised, which is currently 
lacking. 

In addition to hinting at the inadequacy of commonly 
used numerical schemes to yield a sufficiently accurate ex- 
ploration of the energy landscape of XY models from the 
point of view of the determinant criterion, the presence 
of singular, non-isolated stationary points (even after ex- 
plicitly breaking the global 0(2) symmetry by fixing one 
spin) has another relevant consequence. It implies that 
requirements for the validity of the determinant criterion 
(4) itself, as well as of the other theoretical tools devel- 
oped for the study of phase transitions based on station- 
ary points of the energy landscape, are not met by the 
XY Hamiltonian (8). Indeed, all these tools require that 
stationary points are nonsingular and isolated. This is 
typically assumed to be a 'safe' hypothesis once global 
invariances of the Hamiltonian have been removed, but 
our results show that this is not the case. 

This observation may suggest that the application of 
theoretical tools based on the assumption of isolated, 
nonsingular stationary points is hopeless in the case of 
XY models. This is not necessarily true, because a 
way out consists in adding a generic perturbation to the 
Hamiltonian. We have shown in Sec. VII that the singu- 
lar stationary points can be removed by applying generic 
perturbations like (32) or (33). More precisely, the re- 
moval of all the singular stationary configurations has 
been shown for small lattices by the homotopy continu- 
ation method. For larger lattices we have considered a 
sample of Ising stationary configurations, that would be 
singular in absence of the perturbation and that remain 
stationary also in presence of a perturbation of the form 
(33), and we have shown that they become nonsingular 
when the Hamiltonian is perturbed. In previous works 
both the fixing of a single degree of freedom and the ap- 
plication of a perturbation, typically like (33) but with a 
homogeneous field h (see e.g. Ref. [16]), have been consid- 
ered and were thought to be equally effective in removing 
singular, non-isolated solutions. Our results show that 
the global 0(2) symmetry is not the only cause of sin- 
gular solutions, and the 'local' strategy is therefore not 
sufficient for destroying them. 

Remarkably, after switching on even a tiny perturba- 
tion, the rescaled determinant immediately takes on val- 
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ues in the vicinity of the thermodynamic average, far 
from the singular behavior with D = 0. This result tells 
us that the study of the rescaled Hessian determinant D 
carried out in Sec. IV A directly gives us information on 
the behavior of D for the perturbed Hamiltonian (33) 
in the limit of very small external hclds. Since we can 
now safely assume that the perturbed Hamiltonian has 
only isolated singular points, the results shown in Figs. 
1 and 3 should be a faithful representation of what can 
be learned by standard unbiased numerical techniques as 
those employed in this work. 

The presence of families of non- isolated stationary con- 
figurations with zero Hessian determinant (even after 
breaking the global 0(2) symmetry) has interesting im- 
plications reaching beyond the XY models investigated 
in the present paper. Previously, non-isolated stationary 
configurations had already been found in the mean-field 
XY model, but only at a specific value of the energy 
density [16], and also in the globally coupled Kuramoto 
model with homogeneous frequencies [78] (in this con- 
text a continuous family of singular solutions has been 
termed an 'incoherent manifold'). A numerical check by 
means of the homotopy continuation method gave sim- 
ilar results for a variety of other models, including the 
mean-field spherical p-spin model, particles interacting 
via a Lennard- Jones potential and the generalized Thom- 
son problem, details of which will be reported elsewhere. 
Isolated and non-isolated singular solutions often play 
relevant roles also in field theories (see e.g. Ref. [34] and 
references therein). 
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Appendix A: Singular stationary points of the 
three-dimensional XY model 

Similar considerations as for the two-dimensional case 
in Sec. VI motivate the following construction of singular 
stationary configurations in three dimensions, which for 
illustrational purposes is shown here for a lattice of side 
length L = 8. The scheme consists of four different planar 
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(A2) 



(A3) 



(A4) 



All these planar configurations consist of a pattern sim- 
ilar to the two-dimensional configuration (29), but in A 
and C the pattern is shifted one site away from the di- 
agonal. Moreover, C is obtained from A by rotating all 
spins by 7r, and the same is true for D and B. As before, 
the lattice sites marked with gray ^-arrows can be filled 
with an arbitrary 'Ising-typc'-pattcrn of t and J, arrows. 
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Arranging these planar configurations in the sequence 
ABCDABCD results in a stationary configuration on 



a cubic lattice with vanishing Hessian determinant. The 
scheme works in just the same way for larger lattice sizes 
with side lengths that are multiples of 4. 
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